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The formation of higher-order hierarchical systems in star 
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ABSTRACT 

We simulate open clusters containing up to 182 stars initially in the form of singles, 
binaries and triples. Due to the high interaction rate a large number of stable quadru- 
ples, quintuples, scxtuples, and higher-order hierarchies form during the course of the 
simulations. For our choice of initial conditions, the formation rate of quadruple sys- 
tems after about 2 Myr is roughly constant with time at ~ 0.008 per cluster per Myr. 
The formation rate of quintuple and sextuple systems are about half and one quarter, 
respectively, of the quadruple formation rate, and both rates are also approximately 
constant with time. We present reaction channels and relative probabilities for the for- 
mation of persistent systems containing up to six stars. The reaction networks for the 
formation and destruction of quintuple and sextuple systems can become quite com- 
plicated, although the branching ratios remain largely unchanged during the course 
of the cluster evolution. The total number of quadruples is about a factor of three 
smaller than observed in the solar neighbourhood. 
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1 INTRODUCTION 

A sizable fraction of (and possibly all) stars are formed in 
stellar clusters (Lada & Lada 2003, de Grijd et al. 2003, 
Fall et al. 2005). Investigations of starburst regions in 
other galaxies reveal that star clusters are predominantly 
of low mass. For t he small Magellanic cloud and the galax- 
ies M33 and M51, lLamers etahl 120051 find that the initial 
cluster mass function follows a power law distribution with 
a slope of about -2. The minimum mass in their sample 
(M c i ~ 10 4 Mq ) is mainly a result of observational bias. 
The mean mass of star clusters in the solar n eighbour- 
hood is (Af c i) ~ 1000 M (|Kharchenko et al.ll2005h. s u ggest - 
ing that the mass function found by de Griis et al.l (|2003l ) 
should be extended to even lower masses. The vast majority 
(96 ± 2%) of stars of spectral type O do not belong to a 
kno wn association, and may have formed in small star clus- 
ters (|de Wit et al.l 120051) and subsequently been ejected b y 
dynamical slingshots l|Gualandris et al.ll2004lA~arsethl2 005h 
From this it has been arg ued that the m ajority of stars 
are born in small clusters (|Kroupalll995ah . many of which 
may disperse within a few tens of M yr of their formation 
l|Bastian et al.ll2005l ; [Fall et al.ll2005l h In this case, the pop- 
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ulation of Galactic field stars originates mainly from low- 
mass clusters and, as a consequence, the majority of stellar 
multiples in the Galactic field were also born in relatively 
small star clusters. 

The population of binary stars in the solar neigh- 
bourhood has been studie d extensi v ely f o r the field 
|Duauennov fc Mavorl Il99ll : ICouteaul Il993l. Il995l) an d 
for some nearby star clusters ( Kouwenhoven et alj I 



2005). 



Higher-order multiplicities are much more difficult to find 
than binaries. Still, the Multiple Star Catalogue (MSC) con- 
tains 728 systems comprising 3-7 stars. The catalogue is 
claimed to be complet e to a distance of 10 pc from the Sun 
|Tokovininlll997l . Il999h . The majority of the listed systems 
are triples, a nd orbital parameters are p rovided where they 
are available (jSterzik fc Tokovininll2 002') . The MSC contains 
558 triples, 138 quadruples, 25 quintuples, and 7 sextuples. 

We define the multiplicity fraction as the number of 
objects with a given multiplicity divided by the total number 
of single stars and multiples: 

/i = rn/N, 

Here m is the number of objects with multiplicity i and 
N is the sum of the number of single stars (iVs), binaries 
(Nb), triples (Nt), quadruples (ATq d ), quintuples (iVQ n ) and 
sextuples (Afex). Assuming that 15-25% of all systems have 
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multiplicity larger than two (|Tokovininll2004l ), we arrive at 
the following multiplicity fractions: /t = 0.11-0.19, /Qd = 
0.03-0.05, f Qn = 0.005-0.009 and f Sx = 0.001-0.002. (A 
summary of the multiplicities found in the MSC and their 
configuration is presented in Table(2]) 

If most stars are born in relatively small star clusters, 
then the majority of the field population must originate in 
small clusters and the proportions in which single stars, bi- 
naries, triples and higher order multiples occur in the field 
should be reproducible by computer simulations of such star 
clusters. Once these clusters dissolve, the numbers of higher- 
order multiples are no longer affected by cluster dynamical 
evolution, although internal stellar evolution may still affect 
their relative numbers. For example, a triple can evolve into 
a binary when two of its components merge in an unstable 
phase of Roche-lobe overflow. Thus, the relative multiplicity 
fraction is not frozen in when the cluster dissolves. Stellar 
evolution tends to reduce the multiplicity, although these 
effects become less important for older populations. 

Several recent theoretical studies have discussed the for- 
mation of multiple systems. These studies approach the sub- 
ject from two distinct perspectives: (1) stellar dynamical 
models, in which gravitational interactions are computed be- 
tween point-mass stars, and (2) hydrodynamical simulations 
of protocluster evolution. 

Purely dynamical interactions between single stars are 
clearly in effective in producing a suffic iently high fraction of 
binaries (|Aarsetblll97ll ; lAarsethllioql). In addit ion, binary 
orbital periods tend to be too short (|Clarke| [l995). Inclusion 
of hydrodynamical effects can boost the formation rate of 
binaries, but at the cost of reducing further their orbital 
perio ds (|Delgado-Donate et al.l 120031 : iGoodwin fc Kroupal 
2005). Increasing the number of stars in the simulations 
worsens the problem, in the sense that typically only a few 
binaries are formed per star cluster, and high er-order mul- 
tiples are very rare (|Heggie fc Aarsethl Il992h . Simulations 
which start with a large proportion of binaries can repro- 
duce the observed binary frequency l|Kroupall 19951 : lAarsethl 
2004), but studies of the formation of higher-order multi- 
ples in star clusters containing primordial binaries but no 
primordial triples fail to reproduce either the fraction or 
the orbital characteristics o f triples observed in the field 
jPortegies Zwart et al]|2004l ). 

Protocluster evolution, in which gas coagulates to form 
stars or subclusters in the form of hierarchical stellar sys- 
tems, may be able to account for a high proportion of 
binary and possi bly triple systems in young star clusters 
l|Bate et al.ll2003T ). The hydrodynamical breakup of proto- 
stellar cores could be a dominant mechanism for the forma- 
tion of binaries, and possibly triples, in this environment. 
These me thods, however, fail t o produce higher-order mul- 
tiplicities (|White fc G hez 2001). From these studies it would 
seem that star clusters can form from their parent molecu- 
lar cloud with a rich population of binaries and triples, but 
without significant numbers of quadruples, quintuples, etc. 
In that case, these multiples must form during the early 
dynamical evolution of the cluster, by gravitational interac- 
tions between single stars, binaries and triples. 

In this paper we study the formation and reaction rates 
of multiple systems in clusters initially containing ~ 100 
stars, including a sizable fraction of binaries (up to 18%) 
and triples (up to 32%). The simulations are carried out by 



directly integrating Newton's equations of motion to an age 
of approximately 55 Myr, corresponding roughly to the mo- 
ment of dissolution of these clusters. During the evolution, 
cluster members engage in long-lasting hierarchical inter- 
actions involving up to 11 stars. We present the dominant 
interaction channels for multiple stellar encounters in these 
simulations. The small star clusters with primordial binaries 
and triples studied here still under-produce quadruple stellar 
systems compared to the observed fractions of the Galactic 
disc, although the proportions of quadruples, quintuples and 
sextuples relative to one another are roughly consistent with 
observations. 



2 INITIAL CONDITIONS 

The clusters simulated in this study are initialised by select- 
ing positions and veloc i ties fo r 100 objects distributed ac- 
cording to a King ( 1965; 1966) model with Wo = 6. For each 
object we randomly draw a mass between m m i n = O.IMp) 
and m mal = 30Mq from a Kroupa mass function (|Kroupal 
1998), resulting in a total cluster mass of about 45 Mq. Af- 
ter adding binary and triple components, the masses of the 
simulated clusters average about 67 Mq. Finally, we scale 
the velocities of all stars and multiple centers of mass in 
the cluster in such a way that the entire system is in virial 
equil ibrium, in standard AT-body units l|Heggie fc Mathieul 
Il986l ). and scale the cluster to a virial radius r v i T = 0.1 pc, 
corresponding to a half-mass radius of about 0.08 pc. Given 
the rapid expansion of the clusters in the first few Myr 
we consider this a reasonable choice. We call this model S 
(see Tab.[T}. For a ~ 67 Mq star cluster in the solar neigh- 
bourhood, the Jacobi radius in the Galactic tidal field is 
r j ~ 5.5 pc, which is an order of magnitude larger than the 
initial tidal radius of the adopted King model. For clarity 
we therefore neglect the tidal field in our simulations. Note 
that later we describe additional simulations with larger ini- 
tial cluster radii, in which case this assumption may break 
down. 

Starting from model S, we generate model B by ran- 
domly selecting 50 stars, which are converted into binaries 
by adding a companion (secondary) star and orbital param- 
eters. Again, before starting the simulation, the velocities of 
all single stars and binary centres of mass are rescaled so that 
the entire system is again in virial equilibrium. The mass of 
the secondary is randomly chosen from a uniform distribu- 
tion between m m i n and the mass of the primary. The orbital 
parameters are chosen from empirical distributions describ- 
ing approximatel y the inner orbits of hi erarchical triple sys- 
tems in the MSC (|Tokovininll 19971. Il999l ) . The orbital period 
Pbin of the binaries is selected by generating a random vari- 
able X between and 1 and setting 

Pin (X) = 0.09687A" 1 - 0.09677 years. (1) 

The eccentricity for each binary orbit is generated by se- 
lecting another random number Y uniformly in [0,1), and 
defining 

e in (Y) = 1.16 max(0, Y - 0.4). (2) 

If for any binary, the separation at periastron is smaller than 
five t imes the maximum stellar radius f from lEggleton et al.l 
Il989l) . new orbital parameters are selected randomly. 
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Figure 1. Cumulative distribution of the inner and outer or- 
bital period s of 68 observed triple systems (solid curves TTokovinir] 
Il997l . IT999I) . and the distributions (dashed lines) generated from 
Eqs.[TJand[3] 
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Figure 2. Cumulative distribution of the inner and o uter orbital 
eccentricity for 68 observed triple systems (solid curves ; rTokovininl 
1 19071 . 0*9991) . and the distributions (dashed lines) generated from 
Eqs.[2]andH 



These distributions (Eqs.Q] & [2} are presented in Fig- 
ures!]] and [Hand compared with the MSC data. For clarity, 
and due to our lack of understanding of selection effects in 
the discovery of multiple systems, we decided deliberately 
to stay close to the observed distributions for the orbital pe- 
riod and eccentricity. An additional argument for this choice 
is that our simple prescription for the generation of initial 
conditions allows our initial conditions to be easily repro- 
duced. The fact that the binary period distribution does 
not exactly reproduce the observed distribution is therefore 
not a major concern. The other binary parameters (inclina- 
tion, orbital phase an d ascending node) are chosen randomly 
l|Hut fc Bahcalllll983l l. 

Model T is generated by adding a third (outer) star to 
32 randomly selected binaries in model B. The faction of 
triples thus obtained is based on the observed fraction of 
higher order systems in the solar neighborhood. No simula- 
tions were performed with primordial quadruples or higher- 
order hierarchies. The mass of the tertiary star is randomly 
chosen between m m i n and the mass of the inner binary. The 
orbital period Ptrip for the outer orbit of the triple is se- 
lected from a distribution approximating the observed MSC 
distribution of outer orbital periods in triples: 

P out (Z) = 46.8Z' 1 - 45.8 years, (3) 

and the eccentricity for the outer orbit is generated with 

e out (W) = 0.80 W, (4) 

where Z and W are uniformly distributed between and 
1. The other parameters for the outer orbit are chosen ran- 
domly, as for the binary orbits. Final ly, we check the sta- 
bility of the triple using Eq. 90 from Mardling fc Aarsethl 
2001, and new parameters for the outer orbit are selected 
randomly if the triple turns out to be dynamically unstable, 
or if the separation at pericentre between the outer star and 
the binary is less than five times the maximum stellar ra- 
dius. As in model B, the entire system is then restored to 





s 


B 


T 


Ti 


T 3 


n* 


100 


150 


182 


182 


182 


ns 


100 


50 


50 


50 


50 


»B 





50 


18 


18 


18 


TIT 








32 


32 


32 


r v ir [pc] 


0.1 


0.1 


0.1 


1.0 


3.0 


Wo 


6 


6 


6 


6 


6 


N/pc 3 


50k 


50k 


50k 


50 


6 


n Qd 





0.5 


0.78 


0.25 


0.1 


n Qn 





0.1 


0.30 


0.1 


0.01 


™Sx 





0.01 


0.16 


0.05 


0.01 


^runs 


99 


77 


67 


84 


97 



Table 1 . Initial parameters and final fraction of higher-order mul- 
tiples for the various simulations. Each simulation consists of 100 
'particles' which can be single stars, binaries, or triple systems. 
With the adopted number of binaries (hb) and triples (nx) the 
total number of stars n* then ranges between 100 to 182. Each 
simulation was initialised from a King model with Wo = 6 and 
with a virial radius r v ; r = 0.1 pc. Some additional simulations 
were performed with larger virial radii (1 pc and 3 pc). When the 
simulations are terminated at t = 55Myr, each cluster has (on 
average) produced nq d quadruples, riQ n quintuples and ns x sex- 
tuple systems. The final row indicates the total number of simu- 
lations performed for each set of initial conditions, each varying 
only in the original random seed used. 



virial equilibrium by rescaling the velocities of single stars 
and the centre of mass velocities of binaries and triples. 

FigureQ] compares the observed distributions for the in- 
ner and outer orbital period with the distribution generated 
as initial conditions for our simulations. Figure[2] presents a 
similar comparison for eccentricities. Note that our orbital 
separations and eccentricities are uncorrelated with one an- 
other. 
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3 METHODS 

Once the initial realizations are generated, our star clusters 
are evolved u sing the kira integrator in the Starlab soft- 
ware package l|Portegies Zwart et a l. 2001) Q The equations 
of motion are integrated using a fourth-order pre dictor- 
corrector Hermit e scheme (iMakino fc Aarsethlll992l ). using 
block time ste ps (iMcMillanl 1 1986al lbt iMakind 1 199 3 : see also 
lAarsethll2003l for an overview). Stellar and binary evolution 
are modeled with using SeBa |Portegies Zwart fc VerbuntJ 
1 19961 ; iPortegies Zwart et all 1 19971 ). For stable hierarchical 
triples, we follow the internal evolution of the inner bi- 
nary, including Roche-lobe overflow, the effect of stellar- 
wind mass loss, supernovae and the emission of gravitational 
waves. In our prescription within SeBa, the outer star in a 
stable triple cannot initiate a phase of mass transfer to the 
inner binary by Roche-lobe overflow, but the inner binary 
is evolved according to the prescriptions for isolated bina- 
ries. Binaries in higher-order multiples are evolved, so it is 
in principle possible to have a system consisting of six stars, 
in which three binaries each in states of mass transfer, or- 
bit one other. Such a situation, however, is quite rare. All 
simulations were performed on 2.8 GHz workstations Intel 
Pentium 4 processor, and took about 1 hour each. 

In order to clarify some of the results on multiple hier- 
archies we first provide a brief description of our treatment 
and identification procedure for multiple systems. The kira 
integrator employs a dynamic tree structure that represents 
an TV-body system as a mainly 'flat' tree having individual 
stars and the centres of mass of multiple systems as leaves. 
Binary, triple, and more complex multiples are represented 
as binary subtrees below their top-level centre of mass nodes. 
The tree structure determines both how node dynamics is 
implemented and how the long-range gravitational force is 
computed, and also provides a convenient means of identi- 
fying transient structures. 

The tree evolves dynamically according to simple 
heuristic rules: particles that approach 'too close' to one an- 
other are combined into a centre of mass and binary node; 
and when a node becomes 'too large' it is split into its bi- 
nary components. These rules apply at all levels of the tree, 
allowing arbitrarily complex systems to be followed. In prac- 
tice, the term 'too close' is taken to mean that two objects 
approach within roughly the 90° deflection distance for typ- 
ical stellar masses and speeds (see Figure 7.2 on page 421 
of Binney & Tremaine 1987, and also Heggie & Hut 2003, 
where we adopt 9 = tt to compute the 90° deflection dis- 
tance). 'Too large' means that the node's diameter exceeds 
~ 2.5 times this distance. The simulation software prints out 
summary information each time the tree structure changes; 
these data are the basis for the data analysis in this paper. 

Reactions between higher order multiple systems are 
generally chains of events that cause a particular hierarchy 
to change, whether or not the multiple systems are stable. 
The start of a reaction is signalled by a change in the hier- 
archical structure of a multiple. In order to enable a clear 
identification of specific reactions, we introduce a time win- 
dow At, within which a reaction must be completed in order 
to be counted. We set At = 0.01 TV-body time units, which 
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is about the orbital period of a lkT binary and which cor- 
responds to ~550 years in model T. 

A reaction can often be decomposed into a series of 
events in which only a subset of the total number of stars 
participates, each occurring within a smaller time window. 
These events are counted as separate reactions. Reactions 
that last for longer than At are counted multiple times, 
each time the reaction is counted according to the domi- 
nant hierarchy, which sometimes changes continuously. Re- 
actions that last for more than At but that are do not expe- 
rience abrupt changes in hierarchy are counted only once. 
Thus, short-lived resonances are only counted once, but 
those which last for many TV -body time units and change 
configuration may be counted many times. Such long-lived 
resonances, however, are rare. 

In our methodology, we distinguish between multiples 
in terms of the time they remain intact (as just described). 
The survival time of a multiple is taken to be the interval 
between the instant when the multiple forms and the time 
at which components are lost or new ones acquired. If the 
survival time is more than one TV-body time unit (1/2V2 
standard crossing times), we call the multiple 'persistent.' 
Should the survival time be less than this, the multiple is 
called 'transient.' An overview of number of transients and 
persistent multiples created in our simulation model T is 
presented in tables[2] and [3] in § 13.31 

3.1 Evolution of structural parameters 

We concentrate here on the three models S, B and T with an 
initial virial radius r v i r = 0.1 pc, but additional simulations 
have been performed with r v i r = 1 pc and r v i r = 3 pc. The 
latter two models with primordial triples are identified as 
Ti and T3 for r v i r = lpc and r v i r = 3pc, respectively. 
Figure|3]compares the evolution of the bound mass of models 
T, Ti and T3. The means of the Ti and T3 runs are plotted 
as dashed and dotted curves, respectively. The gray shaded 
area in Figure[3] shows the results (mean ±la) for the 67 
simulations of model T. Table[T] summarises the results of 
the various simulations. 

The time at which model T loses half its mass is 
ihm = 15]_7 5 Myr, much shorter than for the other two mod- 
els. A larger initial cluster radius results in a much lower rate 
of mass loss. Mass loss for model T3 is governed by stellar 
evolution. In model T, stellar mass loss is less important, 
and most mass loss is a result of dynamical ejection. The 
sudden divergence between the curves for models Ti and T3 
around 20 Myr is associated with the onset of significant 
dynamical activity in Ti at that time (roughly 2-3 initial 
relaxation times into the evolution). Models S, B, and T 
had comparable overall mass evolution, indicating that an 
initially high proportion of triple systems has little influence 
on the mass loss rate for these clusters; rather, it is mostly 
the relaxation of the cluster that drives dynamical mass loss. 
Simulations that generated more long lived high-order mul- 
tiples (n > 4) did, however, tend to have considerably higher 
mass loss rates. 

The evolution of the cluster half-mass radii for model 
T are presented in Figure [4] Our initial choice of 0.1 pc as 
the initial cluster virial radius seems somewhat small com- 
pared with the o bservations in the open cluster catalogue 
l|Dias et al.ll2002h . However, this is quickly compensated by 
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Figure 3. Evolution of the bound mass for simulations T (grey 
area), Ti (dashed line) and T3 (dots). The grey shaded region 
gives the lcr deviation from the mean for the bound mass of model 
T. For the other two models Ti and T3 we only give the mean 
bound mass; the dispersions are comparable. For model T, the 
upper x-axis shows time measured in JV-body units. 
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Figure 4. Evolution of the half-mass radius for simulations T 
(grey area), Ti (dashed line) and T3 (dots). The grey shaded 
area indicates the extrema (from 67 simulations) measured in our 
simulations of model T. For the other two models Ti and T3, 
we show only the mean bound mass. The dispersions for these 
models are about 0.5 pc for model Ti and 1 pc for model T3. The 
filled circl es show the obser ved star clusters from the open cluster 
catalogue l lDias et al.| [2002). For model T, the upper x-axis shows 
time measured in N— body units. 



the rapid expansion of the cluster, driven by a combination 
of stellar mass loss, dyna mical heating by mu ltiple systems, 
and mass stratification (M erritt et al.l I2004T ). After about 
25Myr the simulated cluster has expanded beyond about 
6pc, which, according to our earlier estimate, would exceed 
the cluster's Jacobi radius in the Galactic tidal field. Al- 
though we ignore the tidal field in our simulations, we argue 
that a cluster which expands beyond this radius should be 
considered dissolved. Some of the smaller observed clusters 
younger than ~ 20 Myr may be satisfactorily reproduced 
by model T. At later times, model T tends to expand too 
rapidly compared to the observed cluster population. How- 
ever, since our models are unlikely to survive this long if the 
Galactic tidal field was taken properly into account the com- 
parison may not be appropriate. We note, however, that the 
expected effect of the Galactic field will be to cause the mass, 
and hence the half-mass radius, of the cluster to begin to 
decline after 25 Myr, possibly improving the agreement be- 
tween the model and the observations, but we do not pursue 
that possibility further here. The initially larger models (Ti 
and T3) may provide somewhat better descriptions of the 
observed clusters at later times. 



3.2 Evolution of multiplicity 

Our main simulations (model T) started with only single 
stars, binaries and triples, but in time a relatively rich popu- 
lation of higher-order multiples formed. During the evolution 
of model T, the binary and triple fractions hardly change, 
while the numbers of stable hierarchical systems consisting 
of 4 stars (nqa), 5 stars (nq n ) and 6 stars (risx) increase 
gradually with time. This is illustrated in Figures[S] and HJ] 
which show the numbers of higher-order multiples (quadru- 
ples, quintuples and sextuples) per cluster as functions of 
time for models T and Ti. 

During the first few million years of model T, the num- 
ber of quadruples increases sharply from zero at birth to 
about 0.35 per cluster at t ~ 2 Myr. After this initial rapid 
increase, the formation rate drops by about a factor of 20, 
and subsequently remains constant until the end of the simu- 
lation. The numbers of quintuples and sextuples show a sim- 
ilar trend of rapid increase at early times and significantly 
slower growth for the rest of the evolution of the cluster, but 
their rates start to increase only after delay times of about 
4 and 10 Myr for the quintuples and sextuples, respectively. 
The numbers of quadruples, quintuples and sextuples at late 
times in model T are quite well approximated by a simple 
least-squares linear fit: 



riQd(t) = 8.0 x 10~ 3 t + 0.34, 
n Qn {t) = 3.4 x 10~ 3 i + 0.11, 
ns*{t) = 1.4 x 10~ 3 i + 0.08. 



(5) 



The formal error on these fits is 4-7%. The formation rate 
drops by about a factor of two for each higher hierarchy 
from quadruples to sextuples. In figure[6] similar trends are 
observable for model Ti. 

It is tempting to interpret the delay in the formation 
of quintuplets and sextuples as the time required to estab- 
lish the channels through which these higher-order systems 
form — we first create quadruples from triples, then quin- 
tuples from quadruples, and finally sextuples from quintu- 
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Figure 5. Numbers of quadruple (solid curve), quintuple (dashes) and sextuple (dots) systems per cluster as functions of time. The 
numbers plotted are valid for the number of systems present at the end of one N— body time unit. These data are averaged over the 
67runs of model T. The dashed lines are fits through the rightmost portion of the data (see Eoj5]l. 
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Figure 6. Numbers of quadruple (solid curve), quintuple (dashes) and sextuple (dots) systems per cluster as functions of time. The 
numbers plotted are valid for the number of systems present at the end of one N— body time unit. These data are averaged over the 
84 runs of model Ti. 



pies, and each stage must wait for a reservoir of lower-order 
systems to form. Unfortunately, as demonstrated in H3.4I 
addition of a single star is not the principal way in which 
quintuples and sextuples form. Rather, the most probable 
formation channels involve triple-binary and triple-triple 
interactions. 

Thus, while it is admittedly difficult to discern the rel- 
evant dynamical details in a single set of runs, we interpret 



Figure[5]as follows. The initial half-mass relaxation time of 
model T is ~ 250 kyr. The relatively rapid initial rise in the 
number of multiples represents the early period of strong dy- 
namical activity when binaries and triples sink to the cluster 
centre and interact, releasing energy and, as a side effect, 
producing some higher-order systems. By 3 Myr the (aver- 
age) cluster has expanded fivefold in radius and its density 
has dropped, reducing the quadruple formation rate mainly 
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by depleting the numbers of single stars. Subsequently, the 
remnant of the cluster is rich in binaries and triples and, 
as the density declines, conditions become more suitable for 
forming and retaining quintuple and sextuple systems. (A 
typical triple in these runs is hard, in the stellar dynamical 
sense, but barely so, making it unlikely that a hierarchical 
system of two such triples could form and survive at the 
initial density of model T.) 

By the end of the simulations (at ~ 55Myr) a total of 
22% (15 out of 67) of the runs have not produced any persis- 
tent hierarchy consisting of 4 or more stars that survived to 
the end of the calculation. A total of 37% (25) of the simula- 
tions have 1 multiple present after 55 Myr, 30% (20) have 2 
multiples, and only 10% (7) of the runs have 3 stable higher- 
order multiple systems present after this time. None of the 
simulations had more than three multiples containing 4 or 
more stars upon termination. Short-lived higher-order sys- 
tems, however, are common in each of the simulations, and a 
total of 7921 higher-order multiples were formed. However, 
most were destroyed before the end of the simulation (see 
§ 13.31 for details, and in particular Tabled- 

Once the cluster disperses, the surviving multiples be- 
come part of the Galactic field, and we may usefully compare 
their numbers with the MSC (§1). Assuming that the field 
was assembled from the remnants of clusters similar to sim- 
ulation T (of which, on average, ~ 100 single stars, binaries, 
and stable multiples remain by the end of the calculation), 
the resultant frequencies of triples, quadruples, quintuples, 
and sextuples are, respectively, 0.25, 0.008, 0.003, and 0.002. 
Thus our particular choice of initial binary and triple param- 
eters under produces quadruples and higher-order multiples 
by a factor of ~ 2-3 compared to the MSC, although the 
relative frequencies of these systems are somewhat encour- 
aging. We now discuss in more detail the breakdown of these 
frequencies within each class. 

Figure|S] gives the the equivalent to Figure[S]but for sim- 
ulation T\. Also in this case the number of multiples in- 
creases quite dramatically in the first few Myr to level off 
at later time, although the overall formation rate of higher 
order multiples is lower than in the models T. Several higher 
order multiples were formed in simulation T3 , but their num- 
ber was rather small and they are omitted from the figure. 

3.3 The hierarchy of multiple systems 

During the simulations, the configurations of the multiple 
systems continually change, as do the numbers of compo- 
nents. In this section we discuss the various hierarchical 
structures, and the frequencies with which they appear in 
the simulations. In § 13.41 we further explore the channels 
through which higher-order multiples are formed and de- 
stroyed. For both this subsection and the next we concen- 
trate on model T. 

In tables[5] and [3] we present an overview of the tran- 
sients and persistent multiples created in simulation model 
T. The multiples which escape the cluster are included in 
this table, as well as the multiples which remain bound. A 
single star in these tables is identified by the letter S. We de- 
note a pair of bound objects by putting parentheses around 
them, with the more massive of the two always to the left. 
For brevity we introduce separate notation for a binary (a 
pair of stars), which we write as B = (S,S). A triple can have 



QU.3ydru.plG conf. 


#MSC 


4L total 


dt 


it P ers - 


dt 


(B,B) 


69 


1171 


5.080 


290 


20.1 


«B,S),S) 


61 


1093 


2.43 


209 


12.1 


((S,B),S) 


1 


945 


0.97 


165 


4.66 


(S,(B,S)) 


•1 


358 


4.45 


41 


38 


(S,(S,B)) 





297 


0.62 


37 


4.0 


total 


138 


3864 




742 




quintuple conf. 


#MSC 


# total 


dt 


# pers. 


dt 




9 


949 


3.08 


202 


13.9 


(R (R 911 


7 


854 


2.16 


145 


11.9 


f(R Rl 91 


5 


142 


2.21 


20 


15 


(KB SI SI SI 


4 


14 


0.33 


2 


2 


(B,(S,B)) 





205 


0.83 


32 


1.5 


<<S,B),B) 





162 


1.83 


29 


9.1 


(S.fB.Bl) 





112 


0.70 


23 


2.8 


C((S,B),S),S) 





12 


0.71 


3 


3 


(S,(S,(B,S))) 





10 


0.43 


2 


2 


((S,(B,S)),S) 





12 


0.13 


1 


1 


(S,((S,B),S)) 





17 


0.07 








(S,((B,S),S)) 





5 


0.01 








((S,(S,B)),S) 





4 


< 0.01 








(S,(S,(S,B))) 





3 


< 0.01 








total 


25 


2501 




460 




sextuple conf. 


#MSC 


# total 


dt 


# pers. 


dt 




3 


88 


6.0 


28 


19 


(fffi Rl 91 Rl 


2 


1 


0.01 








(/R Q"| (R 911 


1 


699 


1.74 


132 


8.53 


( 7R 91 f<Z Rll 


I 


121 


0.74 


11 


6.7 


fY*? Rl (R 911 





381 


0.36 


15 


7.6 







67 


11 


11 


66 


(( (R 91 Rl 91 





16 


0.33 


4 


3 


Cf'9 Rl C9 Rll 





21 


1.2 


1 


30 


to, //R CM R)l 





21 


0.38 


1 


4 


( (R (R 911 91 





17 


13 


3 


71 


(R ( (R 91 911 





14 


0.57 


2 


3 


f(79 Rl Rl 91 





13 


?,.7 








f9 (Y9 Rl Rll 





9 


0.6 


2 


2 


(9 (R (R 9111 





8 


2 


2 


9 


(R (9 (R 9111 





8 


1 


2 


4 


(ffR Rl 91 91 





(5 


0. 1 


1 


2 


^9 (7R Rl 911 





5 


0.5 


1 


2 


(9 (R (9 Rill 

v° t V D t V° i D ) ) ) 





5 


0.3 








(((& S) S) Bl 





4 


0.2 








((B,(S,B)),S) 





4 


0.08 








((S,(B,B)),S) 





2 


0.01 








((S,(B,S)),B) 





3 


0.8 


2 


1 


«S,(S,B)),B) 





1 


0.05 








(B,((S,B),S» 





3 


< 0.01 








(B,(S,(S,B))) 





3 


0.1 








(S,(((S,B),S),S)) 





2 


10 


1 


20 


((((S,B),S),S),S) 





2 


0.09 








(S,(((B,S),S),S)) 





2 


0.01 








(((S,(B,S)),S),S) 





1 


0.8 








total 


7 


1556 




227 





Table 2. Overview of the stellar multiplicities observed in the 
MSC and in our 67 simulations of model T. The first column gives 
this configuration of the multiple. A single star is designated by 
an S, a binary by a B. A binary consists of two single stars and 
could be written as (S, S) = B. The entries are ordered by mass, 
with the more massive component always positioned on the left. 
The parentheses indicate the hierarchy of the multiple system. 
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Table 2 - continued We encounter two stable configurations 
for triple systems, (S, B) and (B, S), 5 stable configurations for 
quadruples, 14 for quintuples, and 29 for sextuples. The second 
column presents the stable configurations listed in the MSC. The 
next two columns give the number of transient and persistent sys- 
tems in our simulation model T, and the mean time (in A^-body 
units) these systems survive before engaging in a new reaction 
(dt). The last two columns give the number of multiples of each 
configuration that survive for at least for one N-body time unit 
(about 55kyr). The final column gives their average lifetimes in 
N— body time units. Note that this average omits systems that 
last for the entire duration of the simulation. The total numbers 
of formations per multiple order (quadruple, quintuple, etc.) are 
presented at the bottom of each column. A summary and a com- 
parison with observed systems is presented in Table|3] 



two configurations: (B,S) if the binary is more massive than 
the star orbiting it, or (S,B) if the outer star is more massive. 
Table[2]draws a distinction between the total number of con- 
figurations and those which are persistent. Table[3] compares 
the simulated hierarchies directly with those observed. 

The most long-lived multiple configuration, on average, 
consists of 4 stars arranged in the classical quartet (B,B), in 
which two binaries orbit one another. A total of 1171 reac- 
tions led to such a configuration. Of these, 290 resulted in 
persisten t syste ms, with an average lifetime of about 20 Myr. 
Aarseth (2004) has commented on the formation of such 
binary-binary systems, and on their longer lifetimes com- 
pared to other, less compact, configurations. Due to their 
smaller cross sections, they are less likely to have a fatal 
encounter with another object than are strictly hierarchical 
systems. The ratios of the numbers of the various quadruple 
configurations to the total number of quadruples in the MSC 
are: 0.50 for (B,B), 0.44 for ((B,S),S), 0.05 for ((S,B),S) and 
(S,(B,S)). The configuration (S,(S,B)), in which the inner 
binary is less massive than the two single stars orbiting it, 
has not been observed. 

To compare these numbers with the results of our sim- 
ulations of model T, we multiply the frequencies at which 
these various configurations are formed by their average life- 
times. The probability of observing any of the above quadru- 
ple configurations during a simulation is then 0.54 for (B,B), 
0.23 for ((B,S),S), 0.14 for ((S,B),S), 0.07 for (S,(B,S)), and 
0.01 for the unobserved (S,(S,B)). These relative ratios are 
comparable to those actually observed. 

In the observations it is often hard to determine the 
masses of the component stars, so we present in Table[3] a 
reduced version of Tabled in which we list the observed 
numbers of hierarchies, but group all systems with the same 
physical hierarchy together. Thus, quadruple systems are 
reduced to just two types (B,B) and ((B,S),S), the latter 
category then contains all possible permutations in which a 
binary is orbited by two single stars: ((B,S),S), ((S,B),S), 
(S,(B,S)) and (S,(S,B)). Among the quintuple configura- 
tions, the strictly hierarchical systems ((B,S),S),S) are vastly 
outnumbered by the triple/binary pairs ((B,S),B) in the 
simulations as well as in the observations. In the observed 
sample, however, permutations of ((B,S),B) are much more 
abundant than in our simulations. The differences between 
the occurrences of observed and simulated multiplicities is 
quite striking. The origin of these discrepancies is not triv- 
ial to understand; it could stem from observational selection 



configuration #MSC # model T 

Common quadruple configurations 
((b,s).s) 69 17.2 ± 0.8 

(b.b) 69 20.2 ± 1.2 



Common quintuple configurations 

(((B,S),S),S) 4 

((B,B),S) 5 1.1 ± 0.2 

<(b,s).b) 9 13.3 ± 0.6 



Common sextuple configurations 
<(b,b),b) 3 4.0 ± 0.6 

(<b,s),(b,s)) 1 4.2 ± 0.4 

«(b,s),s),b) 2 0.9 ± 0.3 



Table 3. Reduced and renormalised version of Table[2] Here all 
possible configurations are summarised without discriminating 
among possible permutations. The number of observed systems is 
given in the second column. The final column shows the number 
of these configurations expected based on the persistent systems 
seen in model T. These numbers are computed on the assumption 
that all stars and multiples in the solar neighbourhood formed in 
clusters similar to model T, and are scaled to the 'survivor' fre- 
quencies presented in the text ( t|3.2ll , assuming a total of 4800 
stars and multiples in the MSC survey. The quoted error is the 
Poissonian uncertainty based on the total number of stable con- 
figurations occurring in simulation T. 

effects or be a result of our choice of initial binary and triple 
parameters. 

3.4 Multiplicity formation and destruction 
reactions 

A gravitational dynamical interaction in which the multi- 
plicity of one object changes can be described in terms of 
creation and destruction reactions. Upon the dissociation of 
a binary, two single stars appear, and when two single stars 
form a bound pair a binary is created. A cascade of such 
reactions then becomes a multiplicity reaction network. It 
is straightforward to extend this description to reactions in 
which more than two stars participate, although the reac- 
tions and the reaction network can become quite compli- 
cated. We now use this perspective to describe the formation 
channels for systems consisting up to six stars. 

The configurations listed in Table[2]and Table|3]are of- 
ten formed in complex interactions involving single stars and 
higher-order systems. In this section we concentrate on the 
individual formation and destruction reactions for multiplic- 
ities of up to 6 stars. 

Figure[7] presents the most common reactions leading to 
the 'creation' (liberation) of single stars due to the gravita- 
tional decay of multiple systems (left), and the correspond- 
ing 'destruction' (consumption) reactions resulting in the 
loss of single stars (right). Each 'ladder' represents a reac- 
tion channel; the symbols S, B and T on the rungs refer 
to single, binary and triple stars, respectively. An arrow is 
drawn from the originating level to the final state of the sys- 
tem. Dashed arrows indicate the reduction of a higher-order 
system to one of lower order — a binary being disrupted or 
a triple being reduced to a single star and a binary. Such 
a reduction could be caused by an encounter with another 
cluster member, by stellar evolution (e.g. a supernova event) 
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Figure 7. Dynamical reactions for the creation (left) and destruction (right) of single stars. In 67 simulations of model T, a total of 
4330 unique creation reactions were identified, versus 3245 unique destruction reactions. Unique meaning that if a particular reaction 
involving a certain set of single star components occurred more than once during a simulation it is counted only once. A reaction results 
in the creation or destruction of at least one single star. In the leftmost creation reaction, for example, two single stars are created from 
the destruction of a single binary via stellar dynamics such as a supernova event of due to the background potential of the star cluster, 
whereas in the second creation reaction one single star and one binary form from a triple. The numbers along the top line show the 
relative frequency in relation to all reactions encountered of this particular reaction. Note that only the most common reactions are 
given, hence these numbers do not add up to unity. The text to the right indicates the state of the system. The following terminology is 
used: S for single stars, B for binaries, T for triples, Qd, Qn and Sx for quadruples, quintuples and sextuples, respectively. 
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Figure 8. Principal reactions leading to the creation and destruction of binaries: 2639 creation and 2628 destruction reactions. 



or, in the case of triple reduction, by the internal dynamics 
of a dynamically unstable system. These different scenarios 
are not represented separately in these diagrams because 
the effect on the multiplicity of the components is the same. 
Solid arrows indicate capture, as when a single star or small 
system combines (typically in the presence of another clus- 
ter member) with another small system to form a higher 
order multiple. Note that, although this presentation might 
suggest that interactions occur in isolation, this is usually 
not the case. Many interactions take place near the core 
of the cluster, and there may be many other stars nearby 
which can carry away small amounts of energy and angular 
momentum. It is sometimes hard to identify the star(s) re- 
sponsible for the binding energy and angular momentum of 
a post-encounter bound pair. Our simulations include stellar 
and binary evolution and allow for stellar collisions to take 
place, but in order not to unnecessarily complicate matters, 
these processes are not displayed separately here. 

The number at the top of each reaction channel indi- 
cates its relative frequency in our simulations. Thus the left 
panel in Figure[7] presents reactions in which single stars are 
created by the ionisation of a binary (40%) or by destruction 
of a hierarchical triple (20%), whereas the right panel de- 
scribes the 'destruction' of single stars by absorption into a 
higher-order system binary (39%) or triple (22%). The sums 
of the creation and destruction frequencies should each be 
unity, but for clarity only the most important reactions are 
shown. The remaining interactions comprise a wealth of low- 
probability, sometimes rather arcane reactions contributing 
to the overall cluster evolution. 

Figuresl8r- 1 121 show analogous reaction channels for the 
formation and destruction of binary through sextuple sys- 
tems. The symbols Qd, Qn and Sx refer to quadruples, quin- 
tuples and sextuples. Note the important role played by 



triples in the formation of higher-order multiples. Quadru- 
ples and quintuples are most commonly formed from a triple 
that absorbs a single star or a binary, whereas sextuples tend 
to form from the combination of two triples. (See also Ta- 
ble[2] where sextuples consisting of two triples orbiting one 
another are relatively common.) We find similar trends in 
the destruction of the higher-order multiples. In particular, 
the decay of sextuples to form two triples is quite striking. 

The reaction frequencies presented in these figures are 
computed for the entire simulation of model T, even though 
one might argue that the early phase of rapid multiple cre- 
ation (t < 2Myr) could produce reactions different from 
those found the later, slower phase of the evolution (Fig. [5}. 
Interestingly, there is no significant difference between the 
types of reactions seen during the early and late phases. 

3.5 The orbital parameters of the hierarchies 

It is not trivial to compare orbital parameters of the ob- 
served multiples with those in the simulations, in part due 
to the large number of parameters and due to severe selec- 
tion effects. Complete orbital solutions are not available for 
any of the multiples listed in the MSC, but for four quadru- 
ples the MSC lists the orbital periods and eccentricities for 
each of the three orbits and also the masses of the four 
stars; those systems are HD 08065+1757, HD 05569+0939, 
HD 11128+3205 and HD 11171-2414. 

Nevertheless we attempt to compare the orbits of the 
multiples in the MSC with those obtained in our simula- 
tions. For clarity we limit ourselves here to the population 
of quadruple systems. Fig. [13] (top panel) shows the orbital 
separations and eccentricities for all quadruples in the MSC. 
The bottom panel shows the same information for persistent 
quadruples at an age of 55 Myr for model T. Small symbols 
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Figure 9. Reactions leading to the creation and destruction of triples: 2700 creation and 2969 destruction reactions. 
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Figure 10. Reactions leading to the creation and destruction of quadruples: 1247 creation and 1229 destruction reactions. 
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Figure 11. Reactions leading to the creation and destruction of quintuples: 724 creation and 760 destruction reactions. 
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Figure 12. Reactions leading to the creation and destruction of scxtuples: 511 creation and 530 destruction reactions. 
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indicate hierarchical systems ((B, S), S), whereas large sym- 
bols show data for the (B, B) configuration. For hierarchical 
quadruples, the outermost, intermediate, and inner orbits 
are represented by small bullets, triangles, and circles, re- 
spectively. For binary-binary systems, the large bullets rep- 
resent the outer orbit, while the large triangles and circles 
indicate, respectively, the inner orbits containing the most 
massive and least massive stars. 

Many of the 'outer' orbits in the observed quadruples 
are missing from the sample in the MSC, which reflects the 
few bullets in Fig. 1131 The outer orbits seem to be cut off at 
around 10 3 AU. The simulated sample of outer orbits, how- 
ever, extends all the way to about 10 5 AU. The few known 
outer orbits in the observed sample of quadruple systems 
strongly suggests that these orbits are simply not observed, 
and that the observational selection effect becomes impor- 
tant at around a few hundred AU. 

There does not seem to be a specific selection effect 
regarding the orbital parameters of the inner binary orbits 
in the 69 (B, B) quadruples listed in the MSC, as for the 
primary as well as the secondary binary a total of 46 are 
observed. The inter-binary orbit, however, is quite heavily 
affected by selection effects, as the MSC lists only 10 binaries 
with known period and eccentricity (the large bullets in the 
top panel of Fig, ll3|l , 

These different selection effects are best seen in the large 
number of missing bullets in the top panel of Fig. 1131 which 
suggests that the outer orbits in hierarchical systems, as well 
as the inter-binary orbits in (B, B) systems, are severely af- 
fected by observational selection effects. The other seem to 
be rather well represented by the observations, with the pos- 
sible exception of the larger number of circularized binaries 
in the observation and some excess of short period eccentric 
binaries in the simulations. This suggests that tidal circu- 
larisation in our simulations is less effective that in nature. 
This should not come as a surprise, as tidal effects are not 
taken into account in higher order systems in which an inner 
binary is perturbed by an outer object, as is often the case 
in these quadruple systems. 

We have performed similar analyses comparing the or- 
bital parameters for the quintuples and sextuple systems 
in the MSC with simulation T, and the results of these are 
quite similar, in the sense that the observational selection ef- 
fects for determining orbital parameters gradually increase 
for larger periods. In these cases we also find that the simu- 
lations tend to overproduce short-period systems with high 
eccentricities relative to the observations, at the expense of 
short-period circular orbits. We decided not to show the 
plots for these higher order multiples, as the quality of the 
observational data declines with increasing multiplicity. 



4 DISCUSSION AND CONCLUSIONS 

We have performed simulations of star clusters which ini- 
tially consisted of single stars, binaries and triples. The ini- 
tial conditions were selected to match the young (< 50 Myr) 
star clusters in the solar neighbourhood, consisting of at 
most a few hundred stars. The initial conditions for the bi- 
naries and triples were selected based on the observed pop- 
ulations, which we transform directly to input parameters. 
For simplicity, observational selection effects are neglected, 
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Figure 13. Semi major axis and eccentricity of quadruples ob- 
served in the MSC (top panel) and those considered persistent 
by the end of simulation T (bottom panel). Each quadruple is 
identified by three sets of orbital parameters; however, for the 
MSC, some parameters are often unknown. The small symbols 
represent hierarchical quadruples ((B, S), S), whereas the large 
symbols are for the (B, B) configurations. For the hierarchical sys- 
tems ((B, S), S), bullets indicate the orbital parameters for the 
outermost orbit, circles represent the inner binary, and triangles 
indicate orbital parameters for the intermediate star orbiting the 
inner binary. For the (B, B) configuration, bullets indicate the 
relative orbital parameters for the two binaries, triangles indicate 
orbital parameters for the binaries containing the most massive 
stars, and triangles represent the less massive binary. Note that 
for the observed sample the orbital parameters are generally not 
all known, whereas for the simulations orbital elements for all 
quadruples are plotted. 



both in determining the distribution functions from which 
we generate our initial conditions, and in our final compar- 
ison with observations. 

After initialisation we run the simulations for about 
55 Myr (corresponding to about 1000 N-body time units, 
or about 350 dynamical crossing times of the system) and 
compare our results with the observed star clusters and mul- 
tiples in the field. Our models with initial half-mass radii 
of about 0.1 pc are consistent with several observed young 
( ^ 20 Myr) star clusters. The rapid expansion of our sim- 
ulated clusters is attributed partly to stellar mass loss and 
partly to dynamical effects. Our simulations fail to explain 
the small ( 2 pc) star clusters at ages £ 20 Myr, although 
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we argue that proper inclusion of tidal effects would mitigate 
this discrepancy. 

The stars in the simulations are initially distributed as 
single stars, binaries and triples. During the evolution, mul- 
tiple systems containing more than 4 stars form as a result 
of strong dynamical interactions among primordial single 
stars, binaries and triples. After the first few million years 
the fraction of higher-order systems remains roughly con- 
stant. 

Our simulation models form significant numbers of hier- 
archical systems consisting of four or more stars when com- 
pared to similar simulations without primordial triples. We 
measure the fractions of systems containing up to 6 stars, 
and find a steady increase in their number. After the start of 
the simulation the number of higher-order multiples rapidly 
increases up to an age of about 2 Myr, after which their net 
formation rate drops by about a factor of 20. The number 
of multiples consisting of 4 stars increases at about twice 
the rate as those containing 5 stars, which again increase at 
about twice the rate of system with 6 stars. 

Among the higher-order multiples that form in our sim- 
ulations, the strictly hierarchical systems are the least likely 
to survive. This is most easily seen in the quadruples, among 
which the configuration consisting of two binaries orbiting 
one another are by far the most common. Among the quin- 
tuples the most likely stable configuration is a triple star 
orbited by a binary. Most sextuple systems consist of two 
triples orbiting one another. 

The relative frequency of stable hierarchies in our sim- 
ulations is generally comparable to those observed in MSC, 
but with some notable exceptions. In four observed systems, 
a massive binary is hierarchically orbited by three single 
stars. In our simulations such sextuple systems are extremely 
rare (relative frequency <C 1%). On the other hand, a sextu- 
ple system consisting of three binaries orbiting one another 
as a hierarchical triple is quite common in our simulations, 
but none are observed in the population of multiple sys- 
tems in the solar neighbourhood. These interesting differ- 
ences may originate from observational selection effects or 
from specific choices in our initial conditions. 

We also present the principal interactions in which mul- 
tiple stars are created and destroyed. Reactions in which 
single stars (~ 50%) or triples (~ 30%) participate are most 
common, simply because such systems are most abundant. 
Most sextuple systems are formed from an interaction be- 
tween two triple systems, and most quadruples form from 
a triple and a single star. The actual relative importance of 
various reactions in the network may be quite different for a 
different choice of initial conditions. The initial fractions of 
stars, binaries and triples may play a crucial role here, and 
a shallower initial density profile may boost the number of 
high order multiples and may even allow for systems con- 
taining larger numbers of stars. We expect, however, that 
multiplicities of 4 or higher will remain relatively rare com- 
pared with systems consisting of fewer stars. 
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